#include <ctime>
#include "TimeIntegratorFactory.h"
using namespace std;

double maxnorm(Vec<double,6> u){
    double M = 0;
    for(int i = 0;i < 3;i++)
        if(abs(u[i]) > M)
            M = abs(u[i]);
    return M;
}

int main(){
    clock_t start,end;
    TimeIntegratorFactory *factory = TimeIntegratorFactory::GetInstance();
    factory->Register("Adams-Bashforth", &AdamsBashforth::Create);
    factory->Register("Adams-Moulton", &AdamsMoulton::Create);
    factory->Register("BDFs", &BDF::Create);
    factory->Register("Runge-Kutta",&RungeKutta::Create);
       
    F f;
    JacobiF ja;
    double T1 = 17.06521656015796;
    double T2 = 19.14045706162071;
    double x1[6] = {0.994,0,0,0,-2.0015851063790825224,0};
    Vec<double,6> u1(x1);
    double x2[6] = {0.87978,0,0,0,-0.3797,0};
    Vec<double,6> u2(x2);
    
    TimeIntegrator *p = NULL;
    string s,filename;
    int order,step,n1,n2;
    Vec<double,6> un,u20;
    double error,e1,rou;
    
    ifstream input("./test/input.data");
    ofstream output1("./run/output_IVP1.data");
    ofstream output2("./run/output_IVP2.data");

    p = factory->CreateObject("Runge-Kutta");
    step = 80000;
    p->T = T2;
    p->k = T2 / step;
    u20 = p->Implement(u2,f,ja,4);
    
    while (!input.eof()){
        input >> s >> order >> n1 >> n2;
        if(input.fail()) break;
        p = factory->CreateObject(s);

        step = n1;
        p->T = T1;
        error = e1 = maxnorm(u1);
        output1 << "IVP1" << "    method:" << s << "    order:" << order << endl;
        output1 << setw(20)<< "time-step" << setw(20) << "solution error" << setw(20) << "convergence rate" << setw(20) << "CPU time" << endl;
        for(int i = 0;i < 6;i++){   
            p->k = T1 / step;
            //filename = "./run/IVP1/" + s + "_" + to_string(order) + "_" + to_string(step) + ".data"; 
            start = clock();
            un = p->Implement(u1,f,ja,order);
            end = clock();
            error = maxnorm(un-u1);
            rou = log(e1/error)/log(2);
            if(i > 0)
                output1 << setw(20)<< step << setw(20) << error << setw(20) << rou << setw(20) << (double)(end - start) / CLOCKS_PER_SEC << "s" << endl;
            else
                output1 << setw(20)<< step << setw(20) << error << setw(20) << " " << setw(20) << (double)(end - start) / CLOCKS_PER_SEC << "s" << endl;
            step = step * 2;e1 = error;
        }
        for(int i = 0;i < 85;i = i + 1)
            output1 << "-";
        output1 << endl;
        
        p->T = T2;
        step = n2;
        error = e1 = maxnorm(u2);
        output2 << "IVP2" << "    method:" << s << "    order:" << order << endl;
        output2 << setw(20)<< "time-step" << setw(20) << "solution error" << setw(20) << "convergence rate" << setw(20) << "CPU time" << endl;
        for(int i = 0;i < 6;i++){
            p->k = T2 / step;
            //filename = "./run/IVP2/" + s + "_" + to_string(order) + "_" + to_string(step) + ".data";
            start = clock();
            un = p->Implement(u2,f,ja,order);
            end = clock();
            error = maxnorm(un-u20);
            rou = log(e1/error)/log(2);
            if(i > 0)
                output2 << setw(20)<< step << setw(20) << error << setw(20) << rou << setw(20) << (double)(end - start) / CLOCKS_PER_SEC << "s" << endl;
            else
                output2 << setw(20)<< step << setw(20) << error << setw(20) << " " << setw(20) << (double)(end - start) / CLOCKS_PER_SEC << "s" << endl;
            step = step * 2;e1 = error;
        }
        for(int i = 0;i < 85;i = i + 1)
            output2 << "-";
        output2 << endl;
        delete p;
    }
    
   
    output1.close();
    output2.close();
    input.close();
    
    step = 24000;
    p = factory->CreateObject("Adams-Bashforth");
    p->T = T1;
    p->k = T1 / step;
    filename = "./run/Euler_24000.data";
    un = p->Implement(u1,f,ja,1,true,filename);
    delete p;

    step = 5000000;
    p = factory->CreateObject("Runge-Kutta");
    p->T = T1;
    p->k = T1 / step;
    filename = "./run/RungeKutta_6000.data";
    un = p->Implement(u1,f,ja,4,true,filename);
    delete p;

}
